#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_statistics.h>

#define C(i) (gsl_vector_get(c,(i)))
#define COV(i,j) (gsl_matrix_get(cov,(i),(j)))

#define BUFSIZE 256

int main(int argc,char *argv[]){
  FILE *input_data=NULL;
  char inp_line[BUFSIZE];
  int n_lines=0,i;
  double inp_freq,inp_flux;
  double *all_freqs=NULL,*all_fluxes=NULL;
  double chisq;
  gsl_matrix *X=NULL,*cov=NULL;
  gsl_vector *y=NULL,*wt=NULL,*c=NULL;
  gsl_multifit_linear_workspace *work=NULL;

  /* load the data */
  input_data=fopen("model_fluxes.txt","r");
  while(fgets(inp_line,BUFSIZE,input_data)!=NULL){
    sscanf(inp_line,"%lf %lf",&inp_freq,&inp_flux);
    n_lines++;
    all_freqs=realloc(all_freqs,n_lines*sizeof(double));
    all_fluxes=realloc(all_fluxes,n_lines*sizeof(double));
    all_freqs[n_lines-1]=inp_freq;
    all_fluxes[n_lines-1]=inp_flux;
  }
  fclose(input_data);

  X=gsl_matrix_alloc(n_lines,4);
  y=gsl_vector_alloc(n_lines);
  wt=gsl_vector_alloc(n_lines);
  c=gsl_vector_alloc(4);
  cov=gsl_matrix_alloc(4,4);
  for (i=0;i<n_lines;i++){
    gsl_matrix_set(X,i,0,1.0);
    gsl_matrix_set(X,i,1,all_freqs[i]);
    gsl_matrix_set(X,i,2,pow(all_freqs[i],2));
    gsl_matrix_set(X,i,3,pow(all_freqs[i],3));
/*     gsl_matrix_set(X,i,4,pow(all_freqs[i],4)); */
    gsl_vector_set(y,i,all_fluxes[i]);
    gsl_vector_set(wt,i,1.0);
  }
  work=gsl_multifit_linear_alloc(n_lines,4);
  gsl_multifit_wlinear(X,wt,y,c,cov,&chisq,work);
  gsl_multifit_linear_free(work);

  printf("1934-638 combined model: log(S) = %g + %g log(f) + %g log(f)^2 + %g log(f)^3\n",/*  + %g log(f)^4\n", */
	 C(0),C(1),C(2),C(3));/* ,C(4)); */

  gsl_matrix_free(X);
  gsl_vector_free(y);
  gsl_vector_free(wt);
  gsl_vector_free(c);
  gsl_matrix_free(cov);

  exit(0);
}
